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Abstract 

Using molecular dynamics (MD) simulations, we find the formation of heaps in a sys- 
tem of granular particles contained in a box with oscillating bottom and fixed sidewalls. 
The simulation includes the effect of static friction, which is found to be crucial in main- 
taining a stable heap. We also find another mechanism for heap formation in systems 
under constant vertical shear. In both systems, heaps are formed due to a net downward 
shear by the sidewalls. We discuss the origin of net downward shear for the vibration 
induced heap. 
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Systems of granular particles (e.g. sand) exhibit many interesting phenomena, such 
as segregation under vibration or shear, density waves in the outflow through hoppers, and 
probably most strikingly, the formation of heap and convection cell under vibration. 1-4 It 
has been known for more than one hundred years that granular particles on the top of a 
vibrating surface will form convection cells and heaps. 5 However, even with many recent 
studies on the subject, 6-11 the exact mechanism for the heap formation is not established. 

Recently, two experimental groups, Evesque et afi and Laroche et al, 7 studied be- 
haviors of granular particles contained in a box, while the whole box is being vertically 
vibrated. They confirm the formation of convection cell and heap. On the other hand, Zik 
et al find convection but no heap. 8 When viewed from above, these boxes are essentially 
squares, making the system fundamentally 3-dimensional. On the other hand, there are 
some studies in 2 dimension (i.e. a line when viewed from above) with fruitful results. 
Using molecular dynamics (MD) simulations of granular particles, Taguchi 9 and Gallas et 
al 10 found convection cells under vibration in 2 dimension. Furthermore, they established 
the fact that the sidewalls are inducing the convection. However, the exact mechanism of 
how the convection is induced by the wall is still not firmly established. Also, they did not 
find any formation of heap. Another breakthrough is the experimental discovery of heap 
formation in 2 dimension by Clement et al. 11 Using monodisperse particles, they found 
that (1) the static friction coefficient must be relatively large in order to induce convection 
and heap, and (2) the heap is formed as particles are being pushed upward by the sidewalls 
(the wall induces convection) along the surface, while there is no significant motion in the 
bulk. The lack of motion in the bulk is probably the consequence of hexagonal packing 
due to monodispersity, and not an essential part of the heap formation. 



The very reason why granular particles can form a stable pile is static friction. More 
precisely, the contact between particles must be able to withstand a finite amount of shear 
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force in order to maintain a pile. We implement static friction in MD simulations of gran- 
ular material using the scheme of Cundall and Strack. 12 In this scheme, one has to apply 
a finite force in order to break a contact between particles. Using the implementations, 
we study heap formations in 2 dimension. First, we present heap formations due to shear 
("the shear induced heaping"), which are intimately connected to "the vibration induced 
heaping." We study the situation that sidewalks are moving vertically in opposite directions 
with constant velocity, thereby creating asymmetrical shear in the cell. Here, the bottom 
plate is not moving. We find the formation of convection and heap. In these simulations, 
the walls are dragging nearby particles, which causes a net flux of particles. This flux is 
inducing convection in the cell, and the convection builds a heap, which is stable due to 
the presence of static friction. We also study the parameter dependence of the formation, 
and find the two static friction coefficients, one between the wall and a particle and the 
other between particles, are the most important. We next study the case that both walls 
are moving down with constant velocity, which causes symmetric shear. We also find a 
convection cell and heaping, whose formations are essentially the same as the asymmetric 
case. Finally, we study the case of vibration induced heaping. We first fix the sidewalls and 
vibrate the bottom plate. We find heap formation and convection for a range of amplitude 
and frequency. Based on several measurements, we propose the following mechanism for 
the formation. The bottom plate is moving up or down during one half of a cycle. The 
density of particles are found to be smaller during the downward phase, which cause the 
shear force by the walls to be larger in absolute magnitude during the upward phase. Over 
one cycle, the net shear force applied by the wall is downward, which cause net downward 
flux of particles near the walls. Therefore, the situation is very similar to the case of the 
symmetric shear. We also study the case of vibrating the sidewalls as well as the bottom 
plate, and find convection but no heap. We discuss possible explanation. 



The force between two particles i and j, in contact with each other, is the following. 
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Let the coordinate of the center of particle % (j) to be Ri (Rj), and f= R\ — Rj. In two 
dimension, we use a new coordinate system defined by the two vectors n (normal) and s 
(shear). Here, n = r/\r\, and s is defined as rotating n clockwise by n/2. The normal 
component of the force acting on particle i by j is 

Fj^i = + a-j ~ \Af /2 - l n m e (v-n), (la) 

where ai (aj) is the radius of particle i (j), and v = df/dt. The first term is the Hertzian 
elastic force, where k n is the elastic constant of the material. And, the constant 7 n of the 
second term is the friction coefficient of a velocity dependent damping term, m e is the 
effective mass, mirrij / {rrii + rrij). The shear component F?^ is given by 

= -l s m e (v-s) - sign(os) min(/c s |os|, ^|Fj^-|), (16) 

where the first term is a velocity dependent damping term similar to that of Eq. (la). The 
second term is to simulate static friction, which requires a finite amount of force (//Fj^J 
to break a contact. 12 Here, fi is the friction coefficient, 5s the total shear displacement 
during a contact, and k s the elastic constant of a virtual spring. There are several studies 
on granular systems using the above interactions. 13 However, only a few of them, 12 ' 14 ' 15 
include static friction. A particle can also interact with a wall. The force on particle i, in 
contact with a wall, is given by Eqs. (1) with cij = oo and m e = rrii. A wall is assumed 
to be rigid, i.e. it is not affected by collisions with particles. Also, the system is under a 
gravitational field g. We do not include the rotation of the particles in present simulation. 
A detailed explanation of the interaction is given elsewhere. 15 

We first consider the situation that systems of granular particles are under constant 
vertical shear. Consider a box of width W and height H. We insert particles at randomly 
chosen positions inside the box, and calculate the trajectories of the particles by a fifth 
order predictor-corrector method. The particles fall by gravity, lose their energy through 
collisions, and fill the box without any significant motion. The parameters we use for the 
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interaction between the particles are k n = 1.0 x 10 6 , k s = 1.0 x 10 4 , 7 n = 1.0 x 10 3 , 7 S = 
and \i pp = 0.2. For the interaction between the particle and the wall, we use k n = 2.0 x 10 6 , 
k s = 1.0 x 10 4 , 7 n = 5.0 x 10 2 , 7 S = 0. The friction coefficient at the sidewall and bottom 
plate are ji pw = 5.0 and 0.2, respectively. The time step is chosen to be 5x 10~ 5 , and gravity 
g is 980. In this letter, CGS units are implied. In order to avoid the hexagonal packing 
formed by particles of the same size, we choose the radius from a gaussian distribution 
with average 0.1 and width 0.02. The density of the particles is chosen to be 0.5. Later, 
we also study the system of monodisperse particles. We then apply a vertical shear by 
pulling the right (left) wall with constant velocity v s = 0.2 (—0.2). In Fig. 1(a), we show 
the system after 80000 iterations of the vertical shear. The slope of the surface of the pile 
increases, and fluctuates around a non-zero value. The mechanism to generate the heap is 
rather simple. Since one pulls the sidewalls with constant velocity, the walls exert shear 
forces to nearby particles. If the force at the wall is sufficiently high, it will induce flow of 
particles in the vertical direction. The upward (downward) flow of particles near the right 
(left) wall, combined with static friction, results in the formation of the heaps. 

We study the effect of parameters on the formation of the heaps. There are quite a 
few parameters in the system. However, most parameters, while their values are chosen 
within reasonable ranges, do not affect the behavior of the system. The key parameters 
are the two friction coefficients n pw and fx pp , and the shear velocity of the sidewalls v s . 
First, we study the effect of n pw . We fix 

A*pp — 0.2, v s — 0.2, and the friction coefficient 
of the bottom plate to be zero. In Fig. 1(b), we show the average angle of the pile (6>) for 
different values of \i pw with W = 3 and the number of particles n = 150. Here, averages 
are taken over time (excluding the transient) and several different runs, where each angle 
is averaged over approximately 5000 points. For small [i pw (0.5 or 1.0), the particles do not 
move significantly during the whole simulations, which results in a zero angle. In order to 
have convective motion and heaping, \x pw should be larger than certain threshold fi pw . The 
existence of a finite threshold \i c pw can be understood as follows. In order to lift particles 

LZ4741 — Jan 5 '93 PAGE 5 



near the right wall, the shear force by the right wall should be larger than the sum of the 
gravitational force and the friction between particles. Since the sum is finite, one needs 
finite n pw in order to maintain the convection. It is still an open question whether the 
transition is the first or the second order, i.e., whether there exists a sudden jump of {9). 

We now fix /i pw = 0.2, v s = 1.0, and study the effect of \i pp . We calculate (9) for 
several values of fx pp , where the averages are taken over approximately 1000 points. Here, 
W = 3 and n = 150. The angle (6) becomes larger for larger values of /j, pp , which may 
results from the fact that the angle of repose is an increasing function of /U pp . 15 We then 
study the effect of v s by fixing \i pw = 5.0, [i pp = 0.2, and change v s . We measure (9) 
for several values of v s between 0.1 and 10.0 with W = 3, n = 150. The measured 
angle is quite insensitive to v s . For example, the angle is 25.8 for v s = 0.1 and 23.6 for 
v s = 10.0. When v s is increased, the pile tries to increase the slope due to larger current of 
particles. On the other hand, increased motion of particles decreases the stablizing effect 
of static friction. These two effects seem to cancel each other resulting to the insensitive 
dependence. Finally, we study the same system using monodisperse particles. The system 
just before we apply the shear is a hexagonal packing of particles with few defects. When 
the shear of v s = 0.2 is applied, the hexagonal packing becomes unstable, and the system 
starts to form a square packing. When the formation of the square packing is completed, 
the particles near the wall can withstand the applied shear with no motion. The square 
packing, which is stable for small values of shear, becomes unstable as v s is increased. 

So far, we have studied the formation of heaps by an asymmetric shear, i.e., the 
sidewalls are moving in the opposite direction. We now consider the case of a symmetric 
shear, where both sidewalls are moving in the same direction. In Fig. 1(c), we show 
the system after 50000 iterations. Here, we use \i pw = 5.0, \i pp = 0.2, and both walls 
are moving down with constant velocity v s = —1.0. The mechanism of generating the 
symmetric heap shown in the figure is essentially the same as that of the asymmetric heap. 
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The shear force induces downward flow of particles near the sidewalls. The flow merges 
together around the center of the cell, and rises to the top of the pile. 

We want to argue that the above "shear induced heaping" is related to the "vibration 
induced heaping." In fact, the above shear geometries are chosen to demonstrate more 
clearly their similarity. We now study the vibration induced heaping, and discuss its 
relation with the shear induced case. We first fix both sidewalls of a box and vibrate the 
bottom plate with amplitude A and frequency /. In Fig. 2, we show the system after 16 
cycles as well as the displacements of the particles over 15 cycles. For this simulation, we 
take W = 10, n = 800, A = 0.190 and / = 20 Hz. The parameters of the interaction of 
the particles and the sidewall are \i = 3.0, fc s = 1.0 x 10 6 . For the interaction between the 
particles, we use ji = 0.5 and k s = 5.0 x 10 4 . All the other parameters are kept the same 
as before. In the figure, one can clearly see a heap and associated convection. 

We now discuss the mechanism for the formation of the heap. In Fig. 2, we show the 
average number of particles c(4>) in contact with one particle for various phases 4> during 
one cycle. Here, 4> is in the unit of 27r, and the averages are taken over 20 cycles. The 
numbers c(4>) are smaller during the downward phase (0.25 < <fi < 0.75) than the upward 
phase. One of the consequences of this "up/down symmetry breaking" is the shear forces 
of the sidewalls are also asymmetric. In Fig. 2, we show the total shear force / s (</>), which 
the right wall applies to the particles, for several values of 4>. The sign of f s is roughly 
opposite to that of the velocity of the bottom. The absolute magnitude of f s is larger 
for the upward phase, and because the particles are more densely packed, the wall can 
exert a larger force. Since the shear force is essentially a drag force for the particles, we 
expect particles move faster vertically during the downward phase, where the shear force 
is smaller. Therefore, there is net downward flux of particles near the sidewalls, which 
results in a convective motion and heaping. In summary, the convection and heaping is 
due to the net current along the sidewalls, which is caused by the net downward shear, 
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which again is a result of "up/down symmetry breaking" of the particle density. 

We consider a similar case when we vibrate the sidewalls as well as the the bottom 
plate. We find that the above mechanism still holds — there is net downward shear and 
flux of particles near the wall, and convection cells. However, we do not find a heap. To 
understand this, we note that the vibration, which is the source of convection and heaping, 
also tends to make an existing heap unstable. The shear force, which is the driving force 
of the convection, is proportional to the relative shear displacement between the walls and 
particles due to the present implementation scheme of static friction. With fixed walls, 
we can obtain rather large net shear force for small A, even when the acceleration of the 
bottom is smaller than g. With moving walls, however, we need the acceleration to be 
at least as large as g to have the relative displacement and net shear. For large A, the 
vibration induces motion of particles along the top surface, which dominates the current 
due to the convective motion. This destroys any existing heap. This observation does not 
necessarily imply that there is no heap formation with moving walls in 2-dimension. It has 
been known that there are some range of A ("windows") for the formation of a heap. 16 
The width of the window is very likely to be dependent on the parameters of the material. 
The window may be very narrow for the range of parameters we are studying. A similar 
parameter problem could also explain why no heap was found in the experiment by Zik 
et al. 8 . There is limitation on the range of parameters we can simulate. For example, k s 
is very large (ideally, infinity) for real material, and larger k s is more effective in creating 
convection. However, we can not simulate values of k s > 1.0 x 10 6 without significantly 
decreasing the time step. The situation is entirely similar with 7 n — we can not increases 
7 n beyond 1 x 10 3 . As a result, it is possible that the heap can be found in these parts of 
the parameter space. 



It has been observed previously that walls are responsible for the convection and/or 
the formation of heaps, and there have been conflicting arguments on the way how the walls 
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induce the convection. We presented here an argument based on measuring properties 
of the system. Our argument is similar to that of Gallas et al 10 in the sense that both are 
based on the shear force that the walls are exerting on the particles. However, the two 
theories have different mechanisms for the generation of the net shear force. 

In conclusion, we find heap formations in two types of systems — one with constant 
vertical shear, the other a with vibrating bottom and fixed walls. Heaps in both systems 
are caused by a net downward shear. Both of the systems are not studied experimentally, 
and there are many interesting quantities to measure. In the vibration induced heaping, it 
would be nice to check for the existence of a net shear by measuring the shear stress of the 
walls. It would be important to study the parameter dependence of the angle of repose. 
In the shear induced case, the further understanding of a parameter dependences of (9) 
(especially, fj, pw ) is necessary experimentally as well as theoretically. Unfortunately, heap 
formation in 3 dimension can not be explained by this mechanism, since it is known that 
heaps can be formed without a boundary in 3-d. The mechanism for 3-d heap formation 
still remains to be understood. 

I thank Hans Herrmann and Michael Leibig for many useful discussions. 
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Figure Captions 

Fig. 1: Shear induced heap formations: (a) Configuration after 80000 iterations of asymmetric 
shear where the right (left) wall is moving up (down) with constant velocity v s = 0.2. 
(b) The dependence of the angle of repose (9) on the friction coefficient of the wall 
H pw . (c) Configuration after 50000 iterations of symmetric shear, where both walls 
are moving down with constant velocity v s = —1.0. 

Fig. 2: Vibration induced heap formation: Configuration after 16 cycles of vibration with 
vibrating bottom plate and fixed sidewalls. Displacements of particles over 15 cycles 
are also shown. In the insets, we also show the average number of contact c((p) and 
the shear force f s (4>) for the different phases </>. 
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